Boundary layer structure in turbulent Rayleigh-Benard convection 
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The structure of the boundary layers in turbulent Rayleigh-Benard convection is studied by means 
of three-dimensional direct numerical simulations. We consider convection in a cylindrical cell at an 
- - - aspect ratio one for Rayleigh numbers of Ra = 3x 10^ and 3x 10^° at fixed Prandtl number Pr = 0.7. 

Similar to the experimental results in the same setup and for the same Prandtl number, the structure 
. of the laminar boundary layers of the velocity and temperature fields is found to deviate from 

■ the prediction of the Prandtl-Blasius-Pohlhausen theory. Deviations decrease when a dynamical 

rescaling of the data with an instantaneously defined boundary layer thickness is performed and 
5— ( , the analysis plane is aligned with the instantaneous direction of the large-scale circulation in the 

closed cell. Our numerical results demonstrate that important assumptions which enter existing 
classical laminar boundary layer theories for forced and natural convection are violated, such as 
the strict two-dimensionality of the dynamics or the steadiness of the fiuid motion. The boundary 
layer dynamics consists of two essential local dynamical building blocks, a plume detachment and a 
post-plume phase. The former is associated with larger variations of the instantaneous thickness of 
, velocity and temperature boundary layer and a fully three-dimensional local fiow. The post-plume 

(~| ■ dynamics is connected with the large-scale circulation in the cell that penetrates the boundary region 

from above. The mean turbulence profiles taken in localized sections of the boundary layer for both 
\ dynamical phases are also compared with solutions of perturbation expansions of the boundary 

' . layer equations of forced or natural convection towards mixed convection. Our analysis of both 

' boundary layers shows that the near-wall dynamics combines elements of forced Blasius-type and 

M-H ' natural convection. 

C/j '. 

Q , PACS numbers: 

•i— I . 

>>. 

r-j . I. INTRODUCTION 

Turbulent Rayleigh-Benard convection can be initiated in a fluid which is confined between a cold isothermal plate 
^-H ' at the top and a hot isothermal plate at the bottom, given a sufficiently strong temperature difference is sustained. 
^ In the turbulent regime, the majority of the heat is carried by convective transport through the layer or cell. Only 
in the vicinity of the top and bottom plates, where the fluid velocities are small, conductive transport takes over 
and becomes important. As in all other wall-bounded flows, boundary layers form. In the present system these are 
boundary layers of the velocity and temperature fields. The structure of these boundary layers turns out to be crucial 
for a deeper understanding of the local and global transport processes as discussed for example in a recent review 
(Ahlers et al. 2009). Furthermore, the boundary layers interact with a so-called large-scale circulation (LSC) that is 
always established in a closed turbulent convection cell. This LSC can take the form of a single roll for aspect ratios 
of order unity or multiple roll patterns for larger ones (du Puits et al. 2007a, van Reeuwijk et al. 2008, Bailon-Cuba 
et al. 2010, Mishra et al. 2011). On the one hand, the LSC is triggered by packets of thermal plumes - fragments 
of the thermal boundary layers which detach randomly from the top and bottom plates into the bulk of the cell. On 
the other hand, the fully established LSC with its complex three-dimensional dynamics can be expected to affect and 
5J] , partly even drive the laminar flow dynamics close to the walls. This interplay has not yet been studied in detail for 
' cylindrical convection cells and provides one central motivation for the present work. 

From a global perspective the heat transport in a turbulent convection cell, which is measured by the dimensionless 
Nusselt number Nu^ is a function of the three dimensionless control parameters in Rayleigh-Benard convection, namely 
the Rayleigh number Ra, the Prandtl number Pr and the aspect ratio P of the convection cell, i.e. Nu = f{Ra, Pr, P). 
Two scaling theories yield different predictions for the turbulent heat transport in convection based on different 
assumptions on the boundary layer structure. While the scaling theory of Shraiman & Siggia (Siggia 1994) is based 
on a turbulent boundary layer with a logarithmic profile for the mean streamwise velocity, Grossmann & Lohse (2000) 
assume laminar boundary layers of Prandtl-Blasius-Pohlhausen type (Prandtl 1905; Blasius 1908; Pohlhausen 1921) in 
order to estimate the boundary layer contributions to the thermal and kinetic energy dissipation rates. Such laminar 
boundary layer evolves in purely forced convection, i.e. for a laminar fiow over a fiat plate. The temperature is treated 
as a passive scalar (Pohlhausen 1921). 

Measuring the boundary layer structure is, however, difficult in laboratory experiments for high-Rayleigh-number 
convection. The reason is that the thickness of the thermal boundary layer, 6t, decreases as the Rayleigh and thus 
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the Nusselt number grow. This thickness is given by 

where H is the height of the convection cell. For a convection flow at Pr ~ the corresponding velocity boundary 

layer will have a similar thickness of 5^ ~ 5t and will thus decrease similarly with increasing Rayleigh number (see 
e.g. Shishkina ct al. 2010). Detailed measurements of boundary lnyci profiles at higher Rayleigh numbers [Ra > 10^) 
require thus large devices such as the Barrel of Ilmenau for the convection in air (du Puits et al. 2007; du Puits et al. 
2010) or high-resolution particle image velocimetry, as possible for convection in water (Sun et al. 2008; Zhou & Xia 
2010). Statistical time-scries analyses of the mean temperature and velocity profiles in the boundary layer yielded 
deviations from the predicted laminar Blasius profiles (du Puits et al. 2007; Zhou & Xia, 2010). A dynamic rescaling 
of the data with respect to an instantaneous boundary layer thickness (which will be explained further below in the 
text) tends to bring it closer to the Blasius prediction in the water experiment by Zhou & Xia (2010). The latter 
result was also confirmed by a series of two-dimensional direct numerical simulations by Zhou et al. (2010, 2011). 
However in both cases, the large-scale circulation is a (quasi-) two-dimensional flow which cannot fluctuate in the 
third direction. 

Du Puits et al. (2007) concluded from their work that the deviations from the Blasius shape arise due to the 
characteristic near-wall coherent structures - so-called thermal plumes which permanently detach from the thermal 
boundary layer. Direct numerical simulations (DNS) by van Reeuwijk et al. (2008a) for Rayleigh numbers up to 10^ 
support systematic deviations from a laminar boundary layer on the basis of an analysis of the friction factor and the 
Reynolds stress budgets. Their DNS showed that the streamwisc pressure gradients have a large magnitude compared 
to Reynolds stresses and are not zero as in the Blasius case. Recall also that the active nature of the temperature 
field is not incorporated in the Prandtl-Blasius-Pohlhausen theory. 

Complementary to the Prandtl-Blasius-Pohlhausen theory for forced convection similarity solutions for natural 
convection are well-known (see e.g. Stewartson 1958; Rotem & Claassen 1969). Here the buoyancy term remains in 
the momentum equation (see below) and is balanced by a wall-normal pressure gradient. The temperature differences 
initiate now fluid motion. Both, purely forced and natural convection, were subject to perturbation expansions 
towards mixed convection which combines forced and natural convection (Sparrow & Minkowycz 1962; Leal 1973). 
This means that either the active role of temperature is inchided as a small-size effect in forced convection or a 
weak outer flow is imposed in natural convection. Hieber (1973) solved numerically the equations which arise from 
perturbative expansions of forced and natural convection. These classical studies are combined with more recent 
efforts to develop two-dimensional boundary layer models for the plume detachment (Fuji 1963; Theerthan & Arakcri 
1998; Puthenveetil & Arakeri 2005; Puthenveetil et al. 2011). The models assume two-dimensional line-like thermal 
plumes with no signiflcant variation perpendicular to the flow plane. 

In this work, we want to resolve the boundary layer structure and its relation to the large-scale circulation for 
Ra > 10^ by means of three-dimensional DNS. We aim at better understanding of the physical reasons for the 
deviations of the boundary layer profiles from the classical Prandtl-Blasius-Pohlhausen and Stewartson theories for 
forced and natural convection, respectively. We, therefore, conduct two long-time DNS of turbulent Rayleigh-Benard 
convection in a cylindrical cell at an aspect ratio F = 1. Step by step it is tested which assumptions of the original 
derivations of the similarity solutions are satisfied or not. Our studies will include analyses of the LSC, the pressure 
gradient fluctuations, the importance of violations of the two- dimensionality of the flow and the active role of the 
temperature at the isothermal walls. The coupling between both boundary layers is also analyzed. We will show that 
actually most of the original assumptions that enter all boundary layer theories are not established in the present 
cellular flow. Furthermore, we relate locally measured turbulence profiles with the results from idealized mixed 
convection boundary layers. 

The outline of the paper is as follows. In the next section, we summarize the numerical model and the equations of 
motion. We present afterwards the boundary layer profiles from the classical time series analysis and the dynamical 
rescaling procedure. The studies are followed by investigations of the large-scale circulation, the pressure fluctuations, 
and time variations of the local boundary layer structure. In section 4 we resolve the dynamics in the boundary layer 
in a small observation window and relate the findings to results of the boundary layer theory of mixed convection. 
We conclude our work with a summary and an outlook. 



FIG. 1: Snapshot of three-dimensional stream hnes in a turbulent convection cell viewed from the top of boundary layer plane. 
The lines are seeded in a horizontal plane inside the thermal boundary layer. Top: Ra = 3 x 10®. Bottom: Ra = 3 x 10^°. 



II. NUMERICAL MODEL 



The three-dimensional Navicr Stokes equations in the Boussinesq approximation are solved in combination with 
an advection-diffusion equation for the temperature field. The system of equations is given by 



+ = -Z + -^ + »9TS,.' (2) 



dui dui dp d^Ui 

dt dxj dxi dx'j 
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where i,j = x,y,z. Here p{x,y,z,t) is the kinematic pressure, Ui{x,y, z,t) the velocity field, T{x,y, z,t) the total 
temperature field, v the kinematic viscosity, and k the diffusivity of the temperature. The dimensionless control 
parameters, the Rayleigh number Ra, the Prandtl number Pr, and the aspect ratio F are defined by 

l^K K H 

Our studies are conducted for F = 1, Pr = 0.7 and Ra = 3 x 10^ and 3 x 10^". Constant a is the thermal 
expansion coefficient, g the gravitational acceleration, AT the outer temperature difference, R the radius and H 
the height of the cylindrical cell. The characteristic length is H, the characteristic velocity is the free-fall velocity 
Uf = \/gaATH . Times are consequently given in units of the free-fall time T/ — H/Uf. The cylindrical geometry 
requires to switch from Cartesian to cylindrical coordinates, (x, y, z) — > (r, 0, z). No-slip boundary conditions for the 
velocity field components , i.e., Ui = 0, hold at all walls. The top and bottom plate are held isothermal at a fixed 
temperatures Tbottom and Ttop, respectively. The side walls are adiabatic with dT/dr = 0. The grid resolutions are 
NrXN^xN:, = 301 X 513 x 360 for i?a = 3 x 10^ and 513 x 1153 x 861 for i?a = 3 x 10^°, where Nr, and 
are the number of grid points in the radial, azimuthal and axial directions respectively. 

The equations are discretized on a staggered grid with a second-order finite difference scheme (Verzicco & Orlandi 
1996; Verzicco & Camussi 2003). The pressure field p is determined by a two-dimensional Poisson solver after applying 
a one-dimensional Fast Fourier Transform (FFT) in the azimuthal direction. The time advancement is done by a 
third-order Runge-Kutta scheme. The grid spacings are non-equidistant in the radial and vertical directions. In the 
vertical direction, the grid spacing is close to Tschebycheff collocation points. The grid resolutions are chosen such 
that the criterion by Grotzbach (1983) is satisfied plane by plane. We define therefore a height-dependent Kolmogorov 
scale as 



1.3/4 

^^(^)= , , ,,1/4 ' (6) 



where the symbol {•) A,t denotes an average over a plane at a fixed height z and an ensemble of statistically independent 
snapshots. Following Emran & Schumacher (2008) and Bailon-Cuba et al. (2010), we define the maximum of the 
geometric mean of the grid spacing at height z by A(z) = maxA^.A^ [\/^^0^r(?')A2(2)]. The thermal boundary layer 
is resolved with 18 grid planes for Ra = 3 x 10^ and with 23 grid planes for Ra = 3 x lO^'^. Thus the recently discussed 
resolution criterion (Shishkina et al., 2010), which would result in 9 and 13 grid planes for the thermal boundary 
layer, is satisfied and over-resolved by almost a factor of 2 in both cases. 

The Nusselt number is found to be Nu ~ 90.32 ± 0.63 for Ra = 3 x 10^ with a standard deviation of 0.7%. The 
second run at i?a = 3 X 10^° resulted in Nu = 189.65 ± 1.5 which gives a standard deviation of 0.8%. The standard 
deviation is determined in the same way as in Bailon-Cuba et al. (2010). We take the Nusselt number plane by plane 
and determine the fiuctuation about the global mean. 

Figure [T] displays instantaneous 3D velocity fields viewed from the edge of the boundary layer close to the bottom 
plate for two Rayleigh numbers. Although a preferential mean flow direction is observable, we see significant deviations 
from a two-dimensionality as visible by the wavy streamlines. With increasing Rayleigh number the streamline plot 
shows more and more textures on an ever finer scale. 



III. BOUNDARY LAYER ANALYSIS 



A. Vertical mean profiles from time series analysis 



Our numerical approach follows the experimental procedure. The latter consists of measuring time-series of the 
three velocity components or temperature at a given point (r, (f), z) in the cell, computing time-averages and repeating 
the measurement for different values of z. The result of such procedures are mean profiles of temperature or velocity. 
In our direct numerical simulation we compute such time series simultaneously for an array of 40 (and 100) points 
starting from z = H. Probe array 1 is at the centerline. Probe arrays 2, 3 and 4 are at r = 0.88i? and — 0,7r/2 
and TT, respectively (see Fig. [ll^a)). We compare the one-dimensional mean profiles for the horizontal velocity V (as 
defined in du Puits et al., 2007) which is given by 



y{r,(t),z,t) = Ju^{r,(j),z,t) -\-ul{r,(l),z,t) , 



(7) 
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FIG. 2: (Color online) Mean profiles of velocity and temperature at Ra = 3 x 10^. (a) Sketch of the four probe arrays with 
measurement locations. Probe array 1 is mounted at at {r,(ji) = (0,0), array 2 at (0.88-R, 0), array 3 at (0.88-R, 7r/2) and array 
4 at (0.88-R, tt). (b) Mean profile of the horizontal velocity V{rj) as defined in ((T]). (c) Mean profile of the vertical velocity 
component Uzijf). (d) Mean profile of the rescaled temperature H'(r;) which is given by ([9]). The solid lines in panels (b) and 
(d) correspond with the classical Blasius and Pohlhausen solutions (Schlichting, 1957). 




FIG. 3: Mean profiles of velocity (a) and temperature (b) at probe array 1 for two different Rayleigh numbers. Quantities are 
the same as in Fig. [2] The solid lines in both panels correspond with the classical Prandtl-Blasius (a) and Pohlhausen solutions 
(b). 



the vertical velocity component Uz and the normalized temperatures s from the top (t) and bottom (b) plates, which 
are defined as 

^ ' ^ ' ^ T{z^ H/2) - Ttop ' ^ ' 

= nr,^^^^^)-n^^-H/'), (9) 
Tbottom — T(z — H/2) 

with the corresponding profiles arising from the Prandtl-Blasius-Pohlhausen theory (see Figs. [2](b)-(d)). Here r] is 
the similarity variable defined in the appendix in The time series contains 57000 data points for Ra = 3 x 10^ 

(and 23000 for Ra — 3x 10^°) at each position of the probe array. This corresponds to 122 (and 58 for Ra = 3 x 10^") 
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FIG. 4: (Color online) Dynamic rescaling of the mean profiles, (a) Dynamically rescaled mean velocity profile for different widths 
of the analysis plane. For comparison the Blasius profile and the profile without rescaling are added, (b) Dynamically rescaled 
mean temperature profile. Again, the Pohlhausen solution and the mean profile without rescaling are shown additionally, (c) 
Comparison of the rescaling of V* as defined in (|f f p with U* (zj) = {ur{r, (pLSC , z,t)\z — ZySv{t))r for a window of width OAR. 
(d) Ratio of velocity magnitudes in the plane, Ur, and perpendicular to the plane, for a window of width OAR. 



free-fall time units Tf. Similar to the laboratory experiments by du Puits et al. (2007) and Zhou & Xia (2010), we 
detect clear deviations from the Blasius and Pohlhausen solutions which are also shown in the figures. Furthermore, 
significant differences among the four profiles can be seen which are caused by the existing large-scale flow in the cell. 
Our profiles at Ra = 3 x 10^ suggest that probe array 4 (and probably array 3 as well) are significantly altered by 
a mean downward motion while probe array 2 is the region of mean upward motion. The mean downward motion 
seems to be connected with an increase of the boundary layer thickness as the data relax much slower to the Blasius 
profile. In section 3.4. we will show that the LSC is on average almost perfectly aligned with the x-axis (0 = 0) for 
the time interval considered in this particular run. In Fig. [31 we compare the data for the two Rayleigh numbers at 
the centerline. The differences between both data sets are very small. 



B. Dynamical rescaling and fluctuations of the boundary layer thickness 

In the next step, we follow the idea of Zhou & Xia (2010) that was applied in their convection experiment in a 
narrow rectangular cell and investigate if a so-called dynamic rescaling of the boundary layer results in mean profiles 
that come closer to the Prandtl-Blasius-Pohlhausen predictions. Similar to the particle image velocimetry in the 
experiment, we analyze the fields in a small planar window. We take this window in the centre of the cylindrical 
cell. Zhou et al. (2011) found that the boundary layer profiles come closer to the Prandtl-Blasius-Pohlhausen case 
downstream the LSC. This plane is in our case additionally aligned for each snapshot with the direction of the 
instantaneous large-scale wind. This direction is determined by the angle 4>lsc which is defined in (|14p . The window 
has a width of 0.02i?, O.IJ? or 0.4i? starting from the centerline of the cell. In order to improve the statistics, we 
conduct this analysis at the top and bottom plates independently for each snapshot. This implies that the large-scale 
flow direction has to be determined separately at both plates. It is known that the large-scale circulation obeys a 
shghtly twisted roll shape (Xi & Xia 2008). 

The instantaneous velocity boundary layer thickness 5^ (t) is determined as the intersection point of the horizontal 
line through the first local maximum of the velocity profile and the tangent to the profile at the plates. The same 
procedure is repeated for the instantaneous thermal boundary layer thickness Sxit). Vertical distances have to be 
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10 20 30 40 50 60 

FIG. 5: (Color online) Fluctuation of the instantaneous thickness of the velocity and thermal boundary layers as formed at 
the center of the top plate of the cell, (top) Ra = 3 x 10^. (bottom) Ra = 3 x 10^°. The dashed lines in both figures mark 
the end of probe array 1 at which the data are taken for both Ra. Thickness values that exceed this height are related with 
velocity profiles that grow gradually from z — and therefore an intersection point with the horizontal line through the first 
local maximum that lies beyond the end of the probe array. 



rescaled with 



Ov(t) OT{t) 

The resulting velocity or temperature profiles follow by (Zhou & Xia, 2010) 

V*{z:) = {V{r,<PLSC,z,t)\z^z:6,(t)),, (11) 
S*(z;) = {E{r,(l)Lsc,z,t)\z ^ z*^5T{t))r ■ (12) 

Here indicates an averaging with respect to r in the plane that is aligned in (j)LSC and the rescaled temperature 
S is taken at the bottom and top, respectively. The corresponding profiles are shown in Figs. HJa) toHl^c). Contrary 
to the experiments by Zhou and Xia (2010) and the two-dimensional DNS by Zhou et al. (2010, 2011) deviations 
to the Prandtl-Blasius-Pohlhausen profiles remain for all window widths and velocities used. A better agreement is 
however observable when the window is chosen to be narrower in radial extension. A further improvement for the 
width 0.4i? is found when the radial component Ur is used instead of V which is defined in Eq. ([7]). For smaller 
windows, however, the agreement with respect to V was better again than for Ur- The deviations for the temperature 
are more persistent which is caused by the plume detachments as we will see in section 4. We verified that the results 
are statistically converged by varying the number of samples. A shift of the window away from the center of the 
plate or a combination of neighboring windows with angles around 4>lsc did not lead to a better agreement with the 
predictions of the Prandtl-Blasius-Pohlhausen theory. The same holds for smaller window sizes than 0.02i?. 

A first significant difference to the previous analysis can, however, be identified immediately. In Fig. lUJd) we 
compare the magnitude of the velocity (ur) in the analysis plane with the velocity (m^) perpendicular to the analysis 
plane with a window of width 0.4_R. It can be seen that the ratio takes a significant non-negligible value in contrast to 
the two-dimensional and quasi- two-dimensional situation, respectively. At z* = 0.5, the height for which the measured 
data start to differ from the theoretical profile, the ratio is grown up to a value of 0.42. This is one important difference 
to the two-dimensional DNS and the quasi-two-dimensional laboratory measurements that gives a first hint of why 
the deviations from the Blasius prediction persist in our geometry. 
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FIG. 6: Cross-correlation ratio as given in Eq. (|13|l for the time series of the instantaneous boundary layer thickness at the 
bottom plate, (a) Ra = 3 x lO" with approximately 57000 data points, (b) Ra = 3 x lO^'' with approximately 23000 data 
points. The outliers from Fig. [S] i.e. the points that exceed the dashed lines, have been excluded from the analysis. 



C. Fluctuating boundary layer thickness 

Figure [5] shovirs time series of both thermal and velocity boundary layer thicknesses obtained from the time series 
at probe array 1. Shorter sequences of same type are obtained from the analysis in the planar observation window 
which is aligned with the instantaneous wind. We can see that both thicknesses fluctuate strongly for both Rayleigh 
numbers. Similar to the two-dimensional DNS of Zhou et al. (2010), the fluctuations of the velocity boundary 
layer thickness are slightly stronger than those of the thermal boundary layer. In particular, we observe rare large 
thickness events for the velocity which can be related to profiles that grow gradually from z = 0. We performed a 
Fourier analysis of both time series and could not detect a characteristic time scale, but a slowly decaying continuous 
spectrum which indicates a chaotic signal. The cross-correlation ratio which is defined as (Zhou et al., 2010a) 

^ mt)~{Sv)t] [ST{t + T)-{ST)t])t 

is plotted in Fig. IHlfor the fluctuating boundary layer thicknesses at the bottom plate. The symbol {■)t denotes a 
time average. Compared to two-dimensional DNS at Pr = 0.7 (Zhou et al., 2010a), the variation of the function g{T) 
is much less regular. In both of our cases the peak is slightly shifted to the left of zero which would indicate that 
variations of the thermal boundary layer cause variations of the velocity boundary layer. The lead time is however 
shorter as the time that we will identify as the time span for a plume detachment. The correlations between both 
fluctuating boundary layers are less pronounced than in the two-dimensional studies. We conclude that such behavior 
is due to the three-dimensional nature of the boundary layer dynamics. We also tried to conduct a similar analysis 
for the data in the small planes which are aligned with the instantaneous large-scale wind. The number of samples 
was, however, too small for a reliable cross-correlation analysis. 

To extract a characteristic time scale from the time series data, we analyze the variations of boundary layer 
thicknesses about their means and to determine the average time intervals oiSy{t) or ST(t) to cross their corresponding 
means. Our data for Ra = 3 x 10^ indicate that this interval for both boundary layers is about Tcross ~ 0.5T/ which 
will turn out to be the time lag for a plume detachment in a local region close to the plate. We also repeated the 
analysis independently for the top boundary layers and reproduced this result. A characteristic variation time of the 
boundary layer is thus 2Tcross ~ Tf. This time can be interpreted as the time at which plumes detach in a local 
region close to the plates (see also our analysis in the next section). It is short when compared to the average loop 
time of Lagrangian tracers in such a cell of approximately 20 Tf (as found in Emran & Schumacher, 2010) which is a 
characteristic loop time of the LSC. 

D. Fluctuations of the large-scale circulation in the convection cell 

Figure [7] displays the direction (or the angle of orientation) and magnitude of the LSC. The orientation the LSC is 
used for the dynamical rescaling of the boundary layer profiles. One can see that the orientation of the mean fiow at 
the same instant is different at the bottom plate compared to the top plate supporting the idea of a twisted circulation 
roU (Funfschilfing & Ahlers, 2004; Xi & Xia, 2008; Xi & Xia, 2008a). The instantaneous direction, 0lsc, and the 
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magnitude of the large-scale circulation, VlsCi are determined by 

<pLSc{to) = (arctan-^i- — ) , (14 

VLScih) = l^\Jux{x, y, zo, toy + Uy{x, y, za, *o)^^^ , (15) 



where the subscript Ar denotes the average over a circular cross-section with r < 0.88i? at zq — St for the bottom or 
zq = H^St for the top plate. Furthermore we show the root-mean-square of the velocity vector which is perpendicular 
to V By. This crossflow velocity vector is determined by the relation vj^ • v — at each point {x, y, zq) S Ar- 

The quantity is given by 



5(^0) = \/ {v\{x,y,za,ta))A, ■ (16) 



It is seen that the circulation is strongly varying in both amplitude and angle. In case of the angle we do observe a fast 
variation of the orientation over a range of approximately 50° . On average the LSC is almost perfectly aligned with the 
a;-axis (0 = 0) along which we have positioned the probe arrays 1,2 and 4. The amount of fluctuations perpendicular 
to the large-scale wind velocity is also significant and reaches up to 50% of Vlsc- The mean magnitude of Vlsc can 
be used to estimate a LSC turnover time by tlsc = Vlsc ^ 2tt{H/2) « 21T/ which is close to the estimate from 
previous Lagrangian studies as mentioned at the end of section 3.3 (Emran & Schumacher 2010). It is also consistent 
with a LSC turnover time of 18 Tf (which corresponds 35 seconds) in the Barrel of Ilmenau. Furthermore, Ahlers et 
al. (2009a) report time scale of 25 seconds from their helium experiment at F = 1/2 that can be converted into 33 
seconds by multiplication with 4/3 for a unit aspect ratio cell. 

On top of the fast oscillation is a very slow drift of the angle (see panels in the upper row). This indicates that a short 
fraction of a very slow precession of the large-scale circulation is monitored. Such a slow-mode can be present since 
the mean orientation of the roll is not locked in one particular direction as being frequently observed in experiments. 
We are however not able to study this slow mode of motion in our DNS since it would exceed our present numerical 
capabilities in terms of the length of the simulation. Better access to this very slow large-scale dynamics would 
require investigations with low-dimensional models (Brown & Ahlers, 2009) or models obtained by proper orthogonal 
decomposition of the turbulence fields (Bailon-Cuba & Schumacher, 2011). 



E. Pressure gradient and temperature fluctuations in the boundary layer 



Figure [5] shows the temperature, the related temperature fluctuations, which are given by 

T'(x,t)=r(x,t)-(T(z))A,t, (17) 
and the magnitude of the horizontal pressure gradient (mid panel in bottom row) , which is given by 

Data are taken at the edge of the boundary layer in the plane at z = 5t at a time instant of the DNS run at 
Ra = 3 X 10^. The contours of 11, which are plotted in units of the logarithm to the base of 10, imply that the 
pressure field is strongly varying in the horizontal plane at this height. More detailed, we display in Fig. [5] the 
quantity 



nc(x,t) = n(x,i)e(n-c), (19) 

with the Heavyside function Q and a threshold C. The pressure field in the incompressible flow limit is directly 
connected with the flow and thus reflects the high spatial (and temporal) variability of the flow including the large- 
scale circulation as analyzed in Fig. [7l 

It has been discussed in detail by Theertan and Arakeri (1998) and Puthenveettil and Arakeri (2005) that the 
horizontal pressure differences are an essential driver of the velocity inside the boundary layer. In Fig. |9]we compare 
vertical profiles taken with respect to time and different horizontal cross sections A — 27rr. Averages of the radial 
component of the pressure gradient {dp / dr) A t and the Reynolds stress {u'^u'^} A,t are shown as examples. The pressure 
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FIG. 7: (Color online) Direction, magnitude and root mean square of the velocity perpendicular to the instantaneous large-scale 
circulation. Quantities are denoted by <j)LSC, Vlsc and V±^rms- Left column is for the top plate and right column for the 
bottom plate. Data correspond to the analysis in Fig. |4]and are determined ai z — St for the bottom plate and sd, z — H — 5t 
for the top plate, respectively. Data are for Ra = 3 x 10^ . The solid horizontal lines indicate the means of the time series. 



gradient component is non- negligible in the boundary layer. Similar to van Reeuwijk et al. (2008a), we compare here 
the ratio 



J^^{dp/dr{z))A.tdz 



[u'j,u'^) At 



(20) 



Note that both terms contribute to the friction factor. Values of 7 = 1.16, 1.77 and 5.21 were obtained for cross 
sections A with radius R, R/2 and i?/5. We thus confirm their finding that this ratio is significant even in a central 
region where the data come closest to the Blasius profiles. We recall that the pressure gradient would be zero in the 
Blasisus case. 

When the spatial support of He is compared with the temperature distribution in the same horizontal plane (see 
bottom row of Fig. [SJ we observe that maxima of 11 are found mostly in the low-temperature voids in between the 
skeleton of plumes, i.e. in regions which are given by (see left panel in bottom row) by T' < or 



T;(x,0=T'(x,Oe(-T'). 



(21) 



Again, we use the Heavyside function Q. In regions of high pressure gradient the horizontal flow will be accelerated 
and piles up local plumes that eventually detach from the boundary layer. The spatial correlation becomes directly 
visible when both thresholded fields lie and T'^ are multiplied as shown in the right panel of the bottom row of Fig. 
m The area covered by these correlated regions is about 11 per cent of the total area and remained nearly constant 
in time, which we verified by a pressure field snapshot analysis over a few free fall time units in case of Ra = 3 x 10^. 

It is also observed from the top panel of Fig. |8l that the plumes are line-like, however with significant thickness 
modulations along their stems. At Prandtl number 0.7 and for the present Rayleigh numbers, diffusion still plays an 
important role in the plume formation. This will to our view also result in a limited applicability of two-dimensional 
plume models in which spatial variations in the third direction are assumed to be small (e.g. Puthenveetil & Arakeri 
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FIG. 8: (Color online) Spatial correlation between the horizontal pressure gradient and temperature. The top figure shows a 
horizontal cross section of temperature T. The three contour plots below show the thresholded temperature fluctuations 
(left figure), the pressure gradient magnitude lie (mid panel), and the product of both (right panel). Data are for Ra = 3 x 10^ 
and taken at 2 = 5t- Pressure gradient magnitude and product are shown in logarithmic units. 




0.4 0.8 1.2 1.6 2 2.4 0.01 0.02 0.03 0.04 



FIG. 9: Vertical mean profiles of the Reynolds stress (a) and the radial derivative of pressure (b) obtained over different 
horizontal cross sections. The full circular cross section A of the cylindrical cell with r = 0.5 (or r = _R) is compared with 
smaller cross sections A that have half and one fifth of the radius, all of which concentric with respect to the center line. Data 
are for Ra = 3 x 10® and in units of U^. 



2005). Similar temperature patterns have been found in Zhou et al. (2007), Shishkina & Wagner (2008), Zhou & 
Xia (2010b) and Puthenveetil et al. (2011) where length, width and aspect ratio of the filaments in this skeleton of 
plumes have been quantified in detail. 

Figure [TU] displays the probability density function (PDF) of the two components of the pressure gradient in two 
planes parallel to the bottom plate. This figure underlines the findings from Fig. |8l The fluctuations of the pressure 
gradient across the boundary layer are highly intermittent as shown by the stretched exponential probability density 
functions of the both components. In Emran & Schumacher (2008), the statistics of the temperature field and its 
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FIG. 10: Probability density functions of the two horizontal components of the pressure gradient. Figures (a) and (c) display 
dp/dr in the center plane of the cell and at the top of the boundary layer. Figures (b) and (d) show r~^dp/d(j>- Data are for 
Ra^3x 10^ 



gradients has been studied in detail. The spatial variations of the temperature as quantified by the statistics of the 
temperature gradient components as well as the thermal dissipation rate were found to obey the strongest spatial 
intermittency in the boundary layer. The intermittency of the pressure gradient field shows qualitatively the same 
behavior, it is enhanced in the boundary layer. 

We summarize our boundary analysis at this point. The numerical data demonstrate that significant differences 
from the classical Prandtl-Blasius-Pohlhausen theory are present in comparison to the two-dimensional case and 
the quasi-two-dimensional experiments. The near-wall flow and temperature structures are three-dimensional and 
unsteady as the large-scale circulation to which the boundary layer dynamics is coupled. This is in line with a 
fluctuating large-scale circulation and the horizontal pressure gradient in the cylindrical cell. 



IV. COMPARISON WITH LAMINAR BOUNDARY LAYERS OF MIXED CONVECTION 



A. Two-dimensional boundary layer theory of mixed convection 



As already discussed in the introduction, the boundary layer in turbulent convection can be considered as mixed 
type, i.e. driven by the natural convection and additionally by the LSC. In the classical boundary layer theory both 
limiting cases have been studied to some extension. These are the purely forced convective flow also known as the 
classical Prandtl-Blasius-Pohlhausen case (Blasius 1908; Pohlhausen 1921) and the purely natural convective flow 
(Stewartson 1958; Rotem & Claassen 1969). For mixed convection, the Boussinesq equations of motion are 
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reduced to the following set of two-dimensional and steady boundary layer equations (Schlichting 1957) 

dux , du^ dp d'^u^ 

VUz^T— = --^ + ^-^-^ ^ (22) 

ox oz ox Oz'^ 

= -^+agT, (23) 
Oz 

dT dT d^T , , 

Ux^ + u^— = K—^, (25) 
Ox oz Oz'^ 

The corresponding dimensionless parameters are the Reynolds and Grashof numbers of the problem which are given 

by 

Rex = ^, Gtx = . (26) 

At the plate (z = 0) the boundary conditions are T = and Ux — = 0. Far away from the plate (z — > oo) it 
follows that T — Too and 

I s f Vao for forced convection 

Ux(Z'^O0) = < n -c ^ 1 ^- (27) 

^ ' 10 for natural convection ^ ' 

In both cases one can define similarity variables r\ and parameters e for the perturbation expansion of mixed convection. 
In agreement with the definitions (l7|) - ([9]) we can proceed as follows. Starting from purely forced convection, the 
expansion follows to (Sparrow & Minkowycz, 1962) 

= fM+^fM^---. (28) 

•'oo 

S(x, .) = ^[!'"]~^°° = B,{^) + e0i(r;) + . . . , (29) 
while starting from purely natural convection, it reads (Stewartson, 1958) 

Ur{x, z) 



Vn{x) 



= %{v) + ^9'i[v) + --- . (30) 



S(x,.) = ^^^'"2^^°° ^Xo(r/)+exi(^) + ..., (31) 

where functions with index represent the unperturbed velocity components or temperature. Furthermore Vn{x) = 
(vg^a^iTyj — TaoYxY^^ . More details are provided in the appendix for completeness. The resulting systems of 
perturbation equations for the boundary value problems can be solved by a shooting method using a 4th-order 
Runge-Kutta scheme (Hieber 1973). 

Figure [11] shows the resulting mean streamwise flow and temperature profiles for the case of Pr = 0.7. The 
perturbation expansion has been carried out to the first order only and curves are plotted for different magnitudes 
of e as given in (1551) . Several aspects can be observed. The boundary layer flow is accelerated if buoyancy effects are 
added to the classical Blasius case as seen in panel (a) of the figure. The same holds if a purely natural convection 
layer is additionally driven by an outer flow such as the large-scale circulation in the present system (seen panel (b) 
of Fig. [TT|) . The imposed outer flow causes a significant variation of the velocity profile. The modifications in the 
temperature are less pronounced. In both cases the resulting mean temperature profiles deviate slightly from the 
unperturbed results. 

Let us now estimate the Rayleigh number dependence of the boundary layer thicknesses in both limiting cases of 
laminar boundary layers. The dependence of the thickness on the Reynolds number in forced convection is given by 
(Schhchting, 1957) 

<5„ - . (32) 

Rex 

By using a scaling relation between Reynolds and Rayleigh numbers for convection at Pr ~ 1 and F = 1, which is 
taken from Ahlers et al. (2009), namely Re ^ RaP '^^ , this results in a Rayleigh number dependence (5^, ~ (5t ~ JSa""'^^ 
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FIG. 11: Vertical profiles of the normalized downstream velocity and temperature resulting from a first-order perturbative 
expansion of forced (panels a, c) and natural convection (panels b, d). The values of the expansion parameter e are indicated 
in the legend. The similarity variables are given by (|34|l . We also show examples for e being unrealistically large in order to 
indicate the deviations better. 



in the purely forced convection case. In a natural convection boundary layer, the Grashof number is substituted with 
the Reynolds number and with the similarity variable rj = zGrJ'^ jx one obtains 

,33, 

Again, if we are interested here in convection with Prandtl numbers around one such that Gr^ ~ -Ra. It follows then 
that (J^ ~ (5t ^ Ra^^'^ which is very close to the forced case. Both scaling estimates suggest that the differences in 
the Rayleigh number dependence of the boundary layer thicknesses are rather small when both limits - natural and 
forced convection - are compared. With only two runs at different Rayleigh numbers at hand, we are not able are to 
conduct scaling laws of the thicknesses with respect to Ra. 



B. Boundary layer dynamics in a small observation window 



The present DNS give us the possibility to zoom into the boundary layer dynamics at higher Rayleigh numbers 
and to test how close the local profiles match with the results of the classical boundary layer theories that we 
just discussed. Out of the comprehensive data record, we have picked two characteristic dynamic sequences of the 
boundary layer structures - a plume detachment event and the post-plume-detachment phase for which the boundary 
layer relaminarizes again. Each of these typical sequences covers a time lag of about 0.45T/ for our data at both 
Rayleigh numbers. We consider them as the two essential building blocks of the boundary layer dynamics. In order 
to make contact to the classical boundary layer theory, we analyze the fields again in a small vertical observation 
plane that is aligned with the instantaneous large-scale circulation. Our observation window has the size of length x 
height equal to 9St x 9(5t for Ra = 3 x 10^ and of 19St x 19St for -Ra = 3 x 10^°. The dense temporal output of the 
data spans 35Tf for Ra = 3 x 10^ and 5Tf for Ra = 3 x 10^° with a time interval of 0.05T/ in both runs. 

A typical plume detachment event is seen in Fig. [T^] where the temperature is shown in the first, the velocity field 
projected into the plane in the second and the out-of-plane velocity component in the third column, respectively. 
The rise of the hot fluid causes strong upward outflow that is connected with the plume detachment. This is in line 
with a strong inflow in the back of the plume due to the incompressibility of the flow. The whole detachment process 
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FIG. 12: (Color online) Sequence of three vertical plane cut snapshots to illustrate the detachment of a plume. The first column 
plots contours of the temperature T, the second is a vector plot of (wr,Wz), and the third plots contours of u^. The pictures 
in the second and third rows lag behind those of the first row by 0.2T/ and OATf, respectively. The color legend (not shown 
here) for the temperature in the first column corresponds to an equidistant color scale between zero (in blue) and one (in red). 
Scalar magnitudes for the azimuthal velocity component are indicated by the color bar (blue for negative and red for positive 
values ) in the right column. The three time instants are number 1, 5 and 9 out of a sequence of nine equidistant snapshots. 
Data are for Ra ^ 3 x 10^ where St/H = 0.0057. 



is accompanied by a cross wind underlining the three-dimensionality of the whole dynamical process. The magnitude 
of the azimuthal velocity is comparable with the amplitudes of V±^rms in Fig. [T] Furthermore, the largest amplitudes 
of the azimuthal velocity component are found to be in line with the largest values of Sv{t) and Sxit). The plume 
detachment is thus one of the dynamical processes that cause the fluctuations of the boundary layer thicknesses. Our 
snapshot analysis also showed that the thickness variations are not significantly delayed with respect to each other 
which is in line with the short lead time for g{T) which we discussed in section ITlI CI The significant azimuthal velocity 
component confirms previous observations by Shishkina & Wagner (2008) that a strong local vorticity vector field is 
aligned with line-like plume ridge. 

The corresponding mean profiles of all velocity components and the temperature are shown in the left column of Fig. 
1131 They are obtained by averaging in the observational window with respect to the radial direction. The detachment 
is accompanied by a deceleration of the radial velocity and strong upward and downward flows into the bulk region 
as already described above. The temperature profiles deviate significantly from the classical laminar boundary layer 
profile (see Fig. [TT]) as the hot fluid parcel leaves the observation area. 

The ambient post-plume-detachment phase is illustrated in Fig. [TH At a first glance the flow and temperature 
fields seem to agree much better with the predictions from the laminar boundary layer theory. However, not too far 
away from the wall still non- negligible upward and downward flows are present. The stratification of the temperature 
field is nearly unperturbed and the azimuthal component is more homogeneously distributed over the window in 
comparison to the detachment phase. This becomes also obvious from the plots in the right column of Fig. [T^ where 
the temperature profiles are much less perturbed than in a plume detachment phase. Nevertheless, even in this phase 
the flow is three-dimensional as we can see from the profiles of the azimuthal velocity component. Recall that the 
observational window in Figs. [T^] and [U has a height of 0.05H. Up to this distance from the wall, the maximum 
magnitude of the mean vertical velocity component is much smaller. 
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FIG. 13: Sequence of nine profiles tiiat correspond with the data of Fig. [12] for the plume detachment phase (left column) and 
with the data of Fig. [14] for the post-plume phase (right column). They are obtained by an averaging in radial (or x-) direction 
over the window that is shown in Fig. 1121 From top to bottom: Radial velocity component Ur/Uf, vertical velocity component 
u^/Uf, azimuthal velocity component u^/Uf and temperature T. 



The radial velocity and temperature profiles in both sequences indicate that the profiles vary strongly, even over 
such a rather short dynamic sequence. The velocity is strongly enhanced in the boundary layer, as it is also resulting 
from the two-dimensional perturbative analysis, such as in the forced case (see Fig. Illf a)). Furthermore, the presented 
data indicate that the large-scale circulation is always strong enough such that the pure natural convection with a 
streamwise velocity that goes to zero, is not established (see Fig. [TlT b)). 

We repeated this analysis for the second run at Ra = 3 x 10^". The qualitative picture remains unchanged for 
both phases, the plume detachment period and the post-plume phase. Note that the mean advection direction of the 
plumes is now opposite. The data are shown in the same way as for the lower Rayleigh number in Figs. [15] and 1171 
As expected, the detaching plumes are more filamented and the boundary layer in the post-plume phase is thinner. 
The amplitude of the azimuthal velocity component remains significant as seen in Fig. 1161 

In both runs the profiles of show the following behaviour in the vicinity of the wall. In the plume detachment 
phase this velocity component changes the sign when moving forward in time from snapshot one to nine. This is 
not the case in post-plume phase. The differences between the temperature profiles for both phases are even more 
pronounced in comparison to the lower Rayleigh number run. 

The time lags of the plume detachment and post-plume phases have been calculated as follows. We take the 
radially averaged temperature field at z ~ 5i5r for each snapshots in the window. If this value exceeds the mean bulk 
temperature it is assigned with a detachment event otherwise it belongs to the post-plume phase. By applying this 
simple procedure, the time series is digitalized. The resulting step function has shorter and longer time periods for 
both dynamical building blocks. The mean time of plume detachment and post-plume is about the same and gives 
about 0.45T/ for Ra = 3 x 10^ where we had a sufficiently long time series. Combining both gives a typical cycle time 
of Tf which is consistent with the 2Tcross from the fluctuating boundary layer thickness in section 3.4. 

In Figs. 1181 and [TOl we try to match the time averaged profiles obtained from the short dynamic sequences with 
the predictions from the mixed convection boundary layer theory including the first-order perturbation. Our profiles 
display again the features that we have detected in the original time series analysis over much longer time intervals 
(see Fig. [3]). However we can now trace the slower increase of the temperature profile clearly back to the plume 
detachment events. Similar connection holds for the velocity profile in post-plume-detachment phase. The local 
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FIG. 14: (Color online) Sequence of three vertical plane cut snapshots to illustrate the phase after the detachment of a plume. 
The first column plots contours of the temperature T, the second is a vector plot of (wr, Uz), and the third plots contours of u^. 
The pictures in the second and third rows lag behind those of the first row by 0.2T/ and 0.4T/, respectively. The color legend 
(not shown here) for the temperature in the first column corresponds to an equidistant color scale between zero (in blue) and 
one (in red) . Scalar magnitudes for the azimuthal velocity component are indicated by the color bar (blue for negative and red 
for positive values ) in the right column. The three time instants are number 1, 5 and 9 out of a sequence of nine equidistant 
snapshots. Data are again for Ra = 3 x 10^. 



dynamical behaviour suggests that the three-dimensional large-scale circulation is now connected to the boundary 
layer section. Inflows from the top of our observation window are observed which cause large variations of the 
velocity profiles. These variations reach the same magnitude as in the plume detachment phase and manifest in the 
deviations for velocity profile {ur)r in the observation plane (see Fig. [T9)) . We have thus shown that the simulation 
data combine elements of forced and natural convections. Neither in the plume detachment nor in the post-plume 
phase the theoretical profiles of both the temperature and velocity fields can be perfectly matched to the data. The 
dynamics close to the walls is always three-dimensional. 

V. SUMMARY AND OUTLOOK 

We have studied the boundary layer dynamics of three-dimensional turbulent Rayleigh-Benard convection in a 
cylindrical cell of aspect ratio one for Rayleigh numbers larger than 10^. Our studies are focussed on the convection 
in air with a Pr = 0.7. The simulations provide access to the full spatial and temporal information in- and outside 
the thermal and velocity boundary layers. 

The large-scale circulation in the cell is varying in its direction and amplitude significantly providing a time- 
dependent driving of the boundary layer dynamics. The fluctuating LSC is in line with a strongly fiuctuating thickness 
of both boundary layers which can be defined from instantaneous snapshots as suggested by Zhou & Xia (2010). When 
these fluctuations are incorporated into a dynamical rescaling, the matching of the mean profiles to the Prandtl- 
Blasius-Pohlhausen theory improves. However, in the present cylindrical cell, deviations from the classical Prandtl- 
Blasius-Pohlhausen profiles will remain, in particular for the temperature. The profiles do also not fit to the other 
limit case, natural convection. 




FIG. 16: Sequence of nine profiles that correspond with the data of Fig. [15] (left column) and Fig. 1171 (rii 
are given in the same units as discussed in Fig. 1131 



ight column). All data 
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FIG. 17: (Color online) Same sequence of data as in Fig. [14] for Ra = 3 x 10^". 
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FIG. 18: (Color online) Matching of the time-averaged profiles of the plume detachment phase with the predictions from 
boundary layer analysis. The profiles are obtained as time averages from Fig. [13] (left column). Data are compared with purely 
forced and natural convection as well as with a corresponding first-order perturbative expansions, (a) Radial velocity in units 
of Uf for forced convection, (b) Rescaled temperature for forced convection, (c) Radial velocity in units of Uf for natural 
convection, (d) Rescaled temperature for natural convection. Data are the same as in Figs. 1121 and [131 (left column). 
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FIG. 19: (Color online) Matching of the time-averaged profiles of the post-plume-detachment phase with the predictions from 
boundary layer analysis. The profiles are obtained as time averages from Fig. [13] (right column). Panels (a)-(d) are as in Fig. 
1181 Data correspond to Figs. 1131 (right column) and 1141 



In the present DNS we aimed at connecting dynamical behavior in the boundary layer with the observed statistics. 
Our analysis found that the boundary layers follow a three-dimensional dynamics in all dynamical phases. This 
conclusion results from investigations of the pressure, the LSC and local dynamic sequences. Pressure gradient 
components and temperature fluctuate strongly and follow non-Gaussian statistics. A significant flow perpendicular 
to a two-dimensional analysis plane is present during detachment, it is also observed in the post-detachment phase. 
These plumes form a line-like skeleton, but are not found to be parallelly aligned. Their detachment is accompanied 
with a significant variation of the boundary layer thicknesses and a cross-wind (azimuthal velocity) with a significant 
amplitude. 

All these observational outcomes violate the assumptions made in deriving the similarity solutions in the classical 
boundary layer theories. Analyses in a point wise probe array as well as in a observational window support our 
findings. This limits also the applicability of two-dimensional plume models and causes to our point of view the 
deviations from both, the classical Prandtl-Blasius-Pohlhausen and the natural convection cases. It can be expected 
that the dynamics in the boundary layer will become increasingly intermittent when the Rayleigh number grows, 
a point that needs to be investigated further. Such an increasingly intermittent behavior would be typical for a 
transitional boundary layer which is ultimately evolving towards a turbulent one at larger Rayleigh numbers. This 
interpretation would also be in line with the DNS results of the presently highest achievable Rayleigh numbers by 
Stevens et al. (2011). They found that the agreement of a dynamically rescaled thermal boundary layer with the 
Pohlhausen prediction worsens when Ra grows. 

One more point: the previous studies by Puthenveetil et al. (2011) as well as the recent experiments by Zhou and 
Xia (2010) suggest that the velocity boundary layer is much less perturbed when the Prandtl number is increased. 
In this case, the thermal boundary layer thickness becomes much smaller than the thickness of the velocity boundary 
layer. Plumes which detach will have a much narrower stem due to decreased thermal diffusion. We expect therefore 
that the agreement with results from the laminar boundary layer theory will improve. This trend might however be 
compensated by an increasing number of fine-scale textures of the turbulent fields for increasing Rayleigh number. 
Our two streamline plots in Fig. 1 suggest this trend. Further comprehensive numerical and experimental studies are 
thus necessary to answer these questions. 
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VI. APPENDIX: PERTURBATIVE EXPANSION OF THE BOUNDARY LAYER EQUATIONS 



We briefly review here the results reported in Sparrow & Minkowycz (1962) and Stewartson (1959). The similarity 
variable is given by 



and the expansion parameter is given by 



zReJ"^ jx for forced convection 
zGr]J^ /x for natural convection 



Grx/ReV^ for forced convection 
Rcx/Grx^^ for natural convection 



(34) 



(35) 



Since the problem at hand is two-dimensional, one uses the stream function instead of the velocity components which 
automatically satisfies the incompressibility condition (1241) . In the forced convection case the following expansions are 
taken 



tp{x,z) = \/vxVo, 



T(x, z) = Too + [Ty, - To, 



resulting for example in the following expressions for the velocity components 



Uxix,z) = 14c 



.m=0 



Uz{x,z) 



14 



m— 1 



(36) 
(37) 

(38) 
(39) 



where primes denote derivatives with respect to rj. The expansion generates in order e° the classical Prandtl-Blasius- 
Pohlhausen equations 

fo' + lfofo = 0, (40) 
O'o+^foO'o = 0. (41) 
The boundary conditions are /o(0) = /o(0) = 0, 6*0(0) — 1 and /q(cxd) = 1, 6*0(00) = 0. The order reads then 

/r + /o7i + lfof[' - \fU'i - + |/^o = , (42) 

/^o - 00, (43) 
O'i + ^/o^'i - ^r^e, + Pr e'J, = . (44) 

The additional boundary conditions follow to /i(0) = /{(O) = 6*1(0) — and /{(oo) = 0i(oo) = ho{oo) = 0. The last 
two terms of (j42p containing ho and as well as Eq. (1431) arise from the pressure term. In natural convection, the 
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expansions are adapted to 



.m=0 



T{x, z) = Too + [Tw - Too) 



m=0 



The order e° was first discussed by Stewartson (1958) and given by 

1 9 9 

/// If f f 1 1 I 

9o + 5^050 ~ g5o5o " g^o + -^VKq 



0, 



- Xo : 



Xo 



3Pr 



-goxo 



0. 



(45) 
(46) 



(47) 
(48) 
(49) 



The boundary conditions are go{0) = 5o(0) = 0' Xo(0) — 1 and .go(oo) = Xo(oo) = ^0(00) = 0. The perturbative 
expansion to mixed convection with order reads 



1 



k[ 



5 

3Pr 



Pr 2Pr 

Xi + —^gox'i + —g'oxi + —^Xogi 



0, 

Xi 

0, 



(50) 
(51) 
(52) 



with 5i(0) = g'i{0) — Xi(0) — xi(oo) = ki{oo) = and .9'i(oo) — 1. Again, fcp and ki arise from the pressure term. 
Equations (|40I) - PH) and (|T7|) - ([5^ have been solved in order to obtain the results displayed in Fig. [11] 
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